################################################
# CALCULATES PERCENTAGE PROTECTED AREA COVERAGE
################################################
# Unit of analysis: Census Units
# Uses: Shapefiles of protected area coverage intersected with census units
# based on code written by Christoph Nolte
# questions to bowydenbraber@gmail.com


library(xlsReadWrite)
library(foreign)

# Load census unit dataset
cu <- read.dbf("CU00M.dbf")[,c(1,4,2)]
colnames(cu) <- c("cu_id", "area", "pop")

# Define categories and years for which coverage is computed
# ALL - All protected areas (conservation units & indigenous lands)
# ALL_XAPA - All protected areas without APA (Area de Protecao Ambiental)
# UC - All conservation units (strict & sustainable) (Unidades de Conservacao)
# UC_XAPA - All conservation units without APA
# PI - Strictly protected areas (Protecao Integral)
# TI - Indigenous lands (Terras Indigenas)
p.set <- expand.grid(year=c(2000,2010), type=c("ALL", "ALL_XAPA", "UC", "UC_XAPA", "PI", "TI"))

for(i in 1:nrow(p.set)) {
  
  # Set variable names (e.g. "a_ALL_2002")
  p <- p.set[i,]
  p.name <- paste(p$type,"_",p$year,sep="")
  p.areafield <- paste("a_",p.name,sep="")
  p.percfield <- paste("p_",p.name,sep="")
  
  # Read respective shapefile
  ol <- read.dbf(paste("Census_spatial/CU_UC/CU00_",p.name,".dbf",sep=""))
  ol <- ol[,c("CU_ID_50M","F_AREA")]
  
  # Aggregate conservation unit area by census unit
  ol <- aggregate(ol,by=list(ol$CU_ID_50), FUN=sum)[,c(1,3)]
  colnames(ol) <- c("cu_id",p.areafield)
  
  # Merge coverage estimate to census unit
  cu <- merge(cu,ol,by.x="cu_id",by.y="cu_id",all=TRUE)
  # Set NAs to zero
  cu[is.na(cu[,p.areafield]),p.areafield] <- 0
  # Calculate percentage coverage
  cu[,p.percfield] <- cu[,p.areafield]/cu$area
}

# Make varnames shorter (for DBF)
colnames(cu) <- sub("20","", sub("APA","",colnames(cu))) #e.g. "a_ALL_XAPA_2002" > "a_ALL_X_02"

write.dbf(cu, paste("CU00_UC_coverage.dbf"))







